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Abstract. The problem of late time instability in time domain integral equations for electro- 
magnetics is longstanding. While several techniques have been suggested for addressing this problem, 
they either require impractically high degrees of freedom in the basis function or an analytical com- 
putation of matrix elements. The authors recently proposed a method that demonstrates stability 
without requiring either of these. The paper, however, does not present theoretical foundations for 
the choice of, or a rigorous error bounds on the approximation of, the bilinear form. This paper 
complements the authors' previous work by presenting a construction of the bilinear form based on 
the minimization of the energy in the system and a proof for the bounds on the approximation. 
We present results on the bounds developed and few sample scattering results that demonstrate the 
stability of the proposed scheme. 
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1. Introduction and examples. While time domain integral equation (TDIE) 
methods have been used to solve electromagnetic scattering problems for a few decades 
now, instability has been a significant debilitating factor, limiting its practical applica- 
tions. Several methods have been proposed to address this issue mathematically, from 
initial work in acoustics [H [3 H] , to more recent work in 2-D [T^] and 3-D [TTJ [TJ [2] 
electromagnetics. While these methods provide a mechanism to construct solutions 
that are stable for long simulation times [5] , they require high degrees of smoothness 
in the basis function for the proofs to be valid [T]. 

Concurrent with the developments in the mathematical framework, there has 
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also been several numerical attempts to solve the problem, viz. obtain stable solu- 
tions to transient electromagnetic scattering using TDIEs [4, 5 , [6j [14j [15] . While these 
techniques have been quite successful in generating stable solutions for a variety of 
challenging geometries, they rely on analytical computation of the matrix elements. 
However, this precludes the use of either acceleration methodologies or higher or- 
der geometry representations, both of which are essential for application to realistic 
geometries. 

Recently, the authors developed a method [5] that relies on a separable expansion 
of the convolution with the retarded potential Green's function. This scheme approx- 
imates the convolution of the Green's function with the temporal basis function in 
any given closed domain using a summation of smooth polynomial functions and has 
been numerically shown to generate stable solutions for a large class of geometries 
that could previously only be stabilized via the exact integration techniques. While 
this method does not suffer from either of the disadvantages of the schemes developed 
in |14] 115] and provided ample empirical evidence of utility, [9] did not provide a rig- 
orous mechanism for truncating the separable expansion. The present work answers 
this need. 

This paper aims to provide a theoretical footing to the scheme developed in [5] 
ITU] . To this end we will (1) construct a bilinear form for the electric field integral 
equation for time domain electromagnetics whose solution leads to a minimization 
of the energy in the system and (2) provide a rigorous mechanism to truncate the 
polynomial expansion developed in [S] [TO]. While we will provide some numerical 
results that demonstrate the validity of our approximations and the stability of the 
scheme developed, we will defer to [9] for more detailed results on stability. 

The rest of this paper will proceed as follows. In Section [2j we formally state the 
problem following which, Section [3] will construct an energy-based bilinear form of the 
TDEFIE, inspired by |4j El [TJ [5] . The next section will describe the new solution 
scheme developed in [5] and provide the theoretical machinery for choosing bounds 
on the expansion. Finally, Section [5] will provide numerical examples that validate 
this technique. 
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2. Variational Formulation based on the energy . Consider a perfect elec- 
trically conducting (PEC) object that resides in free space. Let denote the surface 
of this object that occupies a domain D_, and let D + = IR 3 /Z)_. Assume that every 
point on the surface is equipped with a unique normal h± (r) pointing into D± such 
that 

(2.1) n = n + = —fi-. 




An electric field parameterized by {E 4 (r, t), H 4 (r, t)} and bandlimited to u> max 
is incident on this body. It is assumed that the fields are quiescent for t < 0. The 
incident field induces unknown currents J(r,r)Vr £ that produce scattered fields 
{E s (r, t), H s (r, t)}. We denote the total electric and magnetic fields in D± as E±(r, t) 
and H±(r, t). The problem then is to solve for the current J(r,i) on the surface of 
the scatterer (fi). To this end, we can setup an equivalent problem by setting the 
fields inside the scatterer (D_) to 

(2.2a) E_(r,iO = -E J (r,i) Vrefl_, 

(2.2b) H_(r,£) = -H>,i) VreD_, 

and the fields outside the scatterer (D + ) to 



(2.2c) 
(2.2d) 



E + (r,() = E s (r,()Vrefl_, 
H + (r,() = H s (r,()VreJ)_, 



4 NAIR, PRAY AND SHANKER 

Then the current J(r, t) on the scatterer surface can be written as 



(2.2e) 



J(r,t) = nx [W(r,t) + W(r,t)]. 



Next, inspired by [51[S1[T7], we derive a variational formulation relating the sum of 
the incident field energy interior and the scattered field energy exterior to D. 



3. An energy inspired variational formulation . The electromagnetic en- 
ergy E(t) at an instant t, M 3 can be written as E{t) = E + {t) + E-(t), where E±(t) is 
defined as 



(3.1) E ± (t) = 



1 



sol dr|E ± (r,i)| 2 + Mo I dr'|H±(r,t)| 

Jd± Jd± 



From Poynting's theorem [TB], we have the rate of change of total energy as 



(3.2) 



Of 



E±(t)= [ drE ± (r,t).(n ± xH ± (r,f)) 
Jn 



Combining (3.2) with (2.1) and (2.2), after some manipulation, yields 



(3.3a) 



d_ 

Of 



E(t) = ^ [E+V) + E-(t)} = J drE s (r, t) ■ J(r, t). 



Thus, the total energy in the system at some time T < oo is given by 



(3.3b) E(T)= I dt~E{t)= [ dt [ drE s (r, t) ■ J(r, 



t), 



since E s (r,i) = OVi < 0. 



The scattered fields E s (r, t) and H ,s (r,t) can be related to the surface current 
J(r, t) via integral operators C and K. as 



(3.4a) E s (r,t) = £o{J(r,t)} 



Mo 

47T 



8 t g o {J(r, t)} - J-V / drQ o {V • J(r, r)} 
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and 



(3.4b) 



n x H s (r, t) = K o {J(r, i)} = nxVx-§o {J(r, t)} 



where Q is the retarded potential operator defined for functions J(r, t) as 



(3.4c) 



Go{3(r,t)} = [ dr' [ dt 

Jn J-oo 



r-r' V ' ; 



Following the development in [?], we note that equation (3.3b) suggests the con- 
struction of a space time variational form as 



(3.5) C dt [ dr j(r, t) ■ [E s (r, t)] = C dt [ dr j(r, t) ■ [-E l (r, t)] 

Jo Jn Jo Jn 



for a test function j(r, t) £ L 2 (Vt x [0, T]). We can write this variational form as 



(3.6) 



£ e [j,J] = dt drj(r,t) • -E*(r,t) 
Jo Jo 



where the bilinear form S e [j, J] can be defined using the operator notation introduced 



in (3.41 



(3.7) B e \j,3}= dt / dr j(r,t). £o{J(r,t)} 



the upper limits on the temporal integral can be restricted to T due to the finite 
velocity of propagation. We can construct a similar bilinear form for the magnetic 



field equation (3.4b) which can be written as 



B h \j,J] = dt drj(r,i). nxff(r,() 
Jo Jn 

(3.8) 



dt / drj(r,i)- (X - /C)o{J(r, t)} 
o Jo 



For some parameter a > 0, we can now construct the variational form for the combined 
field integral equation as 



(3.9) B[j,J]= / dt drj{r,t)- -E'(r,() + anxH'(r,i) 

Jo Jo 
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where the bilinear form B [j, J] can be defined as 



(3.10) B%J]= f dt [ dv j(r,i) • \Co{J(r,t)} + a{l-IC)o{3(r,t)} 

Jo Jo. L 

Finally, we define the notation 

(3.11) <u,v>= f dt [ dr u(r,i) • v(r,t) 

Jo Jn 

for functions u, v <G L 2 (f2 x [0, T]). We can now extend ideas obtained in [51|S1[T7] to 
show coercivity of the bilinear form and stability of the solution scheme. To this end 



we have the following relationship between the bilinear form in (3.9 1 and the energy. 

(3.12) B[3, J] = E(T)+ < J, a (2 - JC) o {J} > . 

Combining this with the compactness of the MFIE (K.) operator [TJ, we have the 
following Carding inequality for the coercivity of the bilinear form B. 

Proposition 3.1. Given the bilinear form, £>[j, J], 3 a constant C > such that 

(3.13) B[3,J]+a< J,(/C)o{J}) > > C ||J|| 2 

for norms in L 2 (Q, x [0, T]). 

The proof of the proposition follows naturally from the bounded nature of the JC 
operator and the positivity of the energy. 

4. Numerical evaluation of the bilinear form . Accurate evaluation of the 



bilinear form in (3.9 1 is important to ensure stability [5]. In practice, a discretization 
subspace Vh, T is chosen as that spanned by a set of space-time basis functions of 
the form S(r)T(i), where S(r) are spatial basis functions (usually the Rao-Wilton- 
Glisson or Thomas-Raviart US] spaces), defined on a triangular tessellation of fl, 
say {f^}. T(t) are usually shifted piecewise Lagrange polynomials with support in 
fit C (0, T). Vh, T is indexed on the average size of the spatial and temporal functions, 
i.e. h « diam(tt s ) and r diam(£l t ). While the demonstration of stability developed 
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in [5j could be extended to electromagnetics p] , that is not the focus of this paper. 
Instead, here we focus on the construction of a rigorous mechanism to evaluate the 
bilinear form above. 

The piecewise nature of the Lagrange polynomials implies that accurate com- 
putation of the bilinear form requires identification of the domains of continuity of 
the integrand HO] • It has been shown that the most accurate way to evaluate the 
integral is to compute the correlation between the source and observation domain 
|15j , identify regions of continuity and evaluate the potential integral on each of these 
regions. However, it has also been shown in [H [51 [T3] that evaluating the integral 
accurately over the source domain and collocation over the observation domain is 
sufficient to lead to stable solutions. This technique is tantamount to to finding arcs 
of intersection between arbitrary triangles (on which the source resides) and spheres 
(centered at the point of observation). While tedious, this can be done for piecewise 
flat triangles [HUTS] . 

In [51 [TU] the authors developed a technique to compute the bilinear form numer- 
ically, without the need to identify domains of piecewise continuity of the integrand. 
This was achieved by approximating the space time convolution of the Green's kernel 
with the basis function using a separable expansion in space and time. In this section, 
we will review the expansion and provide a mechanism to truncate the expansion for 
given error criteria. 

We assume a set of functions of the form j(r, t) = S(r)T(i) defined on a simplicial 
tessellation of the domain such that S(r) = Vr ^ £l s and T(t) = Vt ^ Clf The 
bilinear form then involves convolutions of the form. 



(4.1a) 




which can be rewritten as 



(4.1b) 





|r-r'| 
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Fig. 4.1: Domain of the spatial integrals 

where Ii£i s and IIn t denote indicator functions on f2 s and fl t respectively. Let 
B(ca,q8)(l r |) be the spherical shell in space of radii ca and c/3, and H( a t p)(t) be a 
pulse function in time from t = a to t = /3, such that IIq 3 C B( CQC) g)(|r|) and 
II$7 t C H( a ( g)(f). For a given observation point ro and time to, we have, 

(4.1c) 2o(*o,r ) = 

r> 5 (t -a- 

dr'n a ^(r')B( cct ,^)(|ro-r'|) / dfnn t T(t / )H (aiW (to - t')^-, T 

J-oo l r - r I 



Equation (4. Id) can be re-written using the Stone- Wierstrass's theorem as 



(4.2) l Q {t ,r ) = 

dr'lIo.5(r')B (eQi ^ ) (|ro - r'|) / dt'Un t T(t')M. {a ^{to - t') x 



R3 

E ^r^fci v i r ! r ,| — -Pn (M*. + h) , 



2 |ro-r'| 



Note, the polynomials P n {k\{to — t) + k 2 ) and P n (ki l r °~ r I + fc 2 ) in (4.3) are 
continuous over H( a)| g)(to — t') and B( CQ . C| a)(|ro — r'|) respectively, and specifically, 
over the domains of integration Q s and Qf This permits numerical evaluation of the 
integral without having to resort to identifying domains of continuity. While we have 
demonstrated the technique for the source integral, we note that the same scheme can 
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be extended to both source and testing integrals via computation of the geometric 
correlation function defined in |15) . 



4.1. Truncation of the expansion. The expansion in (4.3) is exact. However 



in practice, this procedure is useful only if the summation can be truncated to some 
N. Thus the aim is to evaluate the error 



e N (r,t)= (4.3) 



Tin, (r)S(r)II nt (t)T(i) *. )t B (cQ , C(3) (|r|)H (a ^ (t) 



r - r' 



N 



n ns (r)S(r)n nt (t)T(t) * ait B (caiC/9) (|r|)H (ai/0 (t) £ & 



1 k^rl 



In order to derive a rigorous bound on the error (4.4), we first define the following 



Definition 4.1. Let the temporal Fourier transform ofT(t) be 



(4.4) 



f(u) = T{T{t)} = / dtT(t)e jut , 



and the spatial Fourier transform of a function S(r) be 



(4.5) 



5(A) = HS(r)} = [ drS(r)e jX ' T . 

JR3 



Using these definitions it is possible to re-write equation (4-4-) as 



(4.6) ejv(w,A) = 



.F{II n ,(r)S(r)IIn t (t)T(t)}T { B (ca , c/3) (|r|)H (Q)/3) (t)- 



S (t -t/-l^l 



r -r' 



(2/ + 1) 



P/(fci 



l r o-r'| 



ro-r' 



Pi(h(t -t') + k 2 ) 



Next, we define spatial and temporal bandwidths of interest A m and U) m . These band- 
widths are typically controlled by the bandwidth of the input signal and the size of the 
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object being analyzed. 



Theorem 4.2. Given the definitions above, 3 constants K\, K 2 and M{max{u) m , A m )) 7 
such that \/N > M, 



(4.7) e N (uj,\\\)<-{z t z s ) N 

z$ 



K\ z t z s 



Ko 1 



4^V 2 (1 - zlzlf 47V2 1 _ Z 2 Z 2 



where z t = and z s = 



Proof. The proof is by construction. The Fourier transform 



iV 



(4.8a) 



(21 



1) Pnfate^+ki) 



1=0 



|r - r'| 



Pn(k 1 (t -t') + k 2 ) 



N 



52 . / 0J 



1=0 



11 \ki) 31 U 



Thus, the error is the residual in the sum given by 



(4.8b) e N (oj,\)= J2 ( 2Z + i : 

l=N+l 



8irc _ h*2 
e 3 k i 



From the definition of the spherical Bessel function, we have 



(4.8c) 



3l( z ) = \/^ J i+(i/2)W- 



Using the monotonicity relationship for cylindrical Bessel functions, for alH + 1 < z. 
(4.8d) J J+( i/2)W< (j-^====Y e(^*) VKz-1, 



and combining with (4.8c I, we have 
(4.8e) 



(0 



STABILITY OF TIME DOMAIN INTEGRAL EQUATIONS 



11 



The quality of this bound is illustrated in Fig. 4.2 Using (4.8e I and (4.8b I, the bound 
on the error is 



-N 



(",|A|)< E 

l=N+l 



(21 + i)J± e -^ e (V' 2 -("WM 2 W' 2 -(|A|c/fci) 2 ) 



(4.8f) 



Y 



ui 



. i+ V f,a -(^) 2 / WM^) / 



Equation (4.8fl provides a tight bound on the error. For ease of representation, we 



can relax the bound by writing 



e N (w,\X\)< 

l=N+l 



(21 + y^e-i^eiVP-^'/W+VP-UWW) 

A V 2/ / 



21 



oo 

* E 

l=N+i L 



(4.8g) < --e 1 ^ E 



=AT+1 



(2/ + 1) 



e 2 ^|A| 
4Z 2 



The series in (4.8g) converges very rapidly and in particular, we have 

'e 2 w|A| s r 



(4.8h) 



e N (u,\M) < ^ E 



i=Af+l 



(2/ + 1) 



47V 2 



which is a standard arithmetico-geometric series with the sum 



Cuj -j^i ( e 2 w|A| 

(4 - 8i) w e 1 



AT 



e 2 w|A| 



27V + 3 



(47V 2 - e 2 w|A|) 2 (47V 2 -e 2 u;|A|) 



Further setting z t = an< ^ Zs = 2k!% > we nave 



2fciJV' 



(4.8j) e *(w,|A|)<^- 



K 2 1 



4AT2 



(1 - z t 2 z 2 ) 2 47V 2 1 - z\zl 



The bound on the error derived above implies that for z t < 1 and z s < 1, the error 
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decreases rapidly with zt and z s . In turn, this also implies that for all values uj > Li m 
and |A| > |A| m , the error rapidly decreases with N, as long as N > 2*max(j^ 1 , j^)- 
Indeed, in practice, it is sufficient to pick N > 2 * min( 5 ^ a , )• We utilize this 
result, to construct an efficient scheme to integrate the bilinear form derived in Section 
[3] In the next section, we will demonstrate numerical results obtained using this 
technique. 

5. Numerical Results. First, we demonstrate the construction of the spatio 
temporal bandwidths, and corresponding truncation limit N. For simplicity and ease 
of demonstration, we will derive the truncation limits in temporal frequency; the 
spatial frequency bounds are identical. 

Consider a band limited signal that is convolved with a Lagrange polynomial of 



order 1. This convolution is then approximated using expression in equation (4.3), 
truncated to different values N. Figure |5.1a| shows the spectrum of the true con- 
volution overlaid on approximations for different N. The domain of the expansion 
M.r a< p\(t) is chosen to be 3 At in width, and the temporal basis functions are chosen to 



be supported in ui t = 2A t . Figure 5.1b shows the convergence in the norm of the error 
within the band as a function of N. As is clear from both the figures, the expansion 
can be efficiently truncated with minimal and controllable loss in accuracy. 

Next, we present results that validate the full solution scheme developed in [3] 
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Fig. 5.1: Numerical realization of truncation bounds 



using the expansions derived in this paper. To this end, we consider scattering from 
a perfectly conducting sphere of radius 1.0m, due to an incident plane wave. The 
incident field is a plane wave of the form 



(5.1) E l (r,i) = ucos(27r/ t)e 



- (t-r-k / c-t p y /2a 2 



where u denotes the polarization vector, fo the center frequency, and k the direction 
of propagation. The values a and t p are calculated as a = 3/(2irB) and t p = 6a, 
where B denotes the bandwidth of E 1 in Hz. The incident power is calculated to 
be approximately 160 dB below the peak at f max = fa + B and f m in = fo — B. 
The sphere is discretized using curvilinear triangular elements (as described in [3]) of 
order g = 3. Temporal basis functions are first order Lagrange polynomials and the 
expansion of the spatio-temporal convolution is truncated using the criteria derived 
in this paper, resulting in a maximum N = 7. The incident field is wide band with a 
center frequency of 150MHz and bandwidth of 149.9MHz. It is seen that the current 



on the surface of the sphere is stable for 20, 000 time steps (Figure 5.2a) with a = 0.5. 
Similarly, the far field observed due to this current, normalized to the incident field, 
can be used as a metric for comparison against an analytical frequency domain result. 
The excellent agreement with the analytical result at three different frequencies across 
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the bandwidth validates the technique (Figure 5.2b I. 






V 



(a) Current on nacelle surface, for early time,(b) Current on nacelle surface for long run time 
comparison against standard MOT scheme showing late time stability of proposed scheme 

Fig. 5.3: Construction of stable solutions on a nacelle geometry 



Finally we present a result that has proven challenging to stabilize using tradi- 
tional MOT based schemes. We compute scattering from a nacelle geometry using 
the proposed scheme. The nacelle is constructed as a cylindrical shell with one end 
closed. The outer radius of the cylinder is 160cm and the inner radius is 100cm. The 
length of the cylinder is 300cm. To be realistic, the edges of the cylinder are filleted 
at a 20° angle. The incident field is chosen to have the same functional dependence 



as in (5.1 1 with /q = 1GHz and B = 0.99GHz. The incident field is polarized along 



9 and directed along 8 = 90° and 4> — 5°. The maximum N chosen across the ge- 
ometry is 7 and the minimum is 2. The simulation is run for 100, 000 time steps, 
corresponding to 753 transits. While the stability of the scheme is proven for a > 0, 
we have chosen a = in this case to illustrate the robustness of the approximation 
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in terms of providing stable solutions for the more challenging case of the electric 
field integral equation. The early time data in Figure |5.3a| shows the current on the 
surface of a nacelle geometry computed using the proposed scheme compared against 
the scheme developed in [14]. The figure not only validates the method but clearly 
indicates the instability of the standard MOT schemes. Late time data in Figure [B".3b| 
demonstrates the capability of the scheme developed in the paper to generate stable 
solutions for extremely long simulation times (753 transits) for problems of scattering 
from challenging geometries. 

6. Conclusions. In this work, we have presented a technique to numerically 
construct a stable solutions to time domain integral equations for electromagnetic 
scattering from perfectly conducting structures. We have derived a bound on the 
truncation error and provided a mechanism for numerically estimating the error. The 
technique has been numerically validated on a canonical structure using higher order 
tessellations (TUJ [TT] and shown to be applicable to generate stable solutions using the 
EFIE for more challenging structures that defy stabilization via traditional methods 
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